Jerry Greenough
Measuring Mesh Quality
This case study demonstrates the measurement of the quality of a large finite element mesh. The measurment is designed to detect quadrilateral elements that have geometric characteristics that deviate from the perfect square. In practice, the majority of finite element meshes have such imperfections. However, it has been shown that if the element shapes deviate beyond certain known bounds, they have the potential to impact the accuracy of stress and deflection predictions generated by a plane-stress finite element analysis. It is therefore wise to measure the quality of each element in a finite element mesh against a number of well-established criteria. For a mesh of quadrilateral elements, the most commonly used criteria are:
The implementation (in C++) of a mesh quality checker is demonstrated here using the first three criteria. A sample mesh (representing a plate with a hole) containing 73344 nodes and 72064 quadrilaeral elements is used to test the implementation. See Figure 2 for a picture of one quadrant of the sample mesh. The mesh is large enough to demonstrate the benefits of employing a parallel multi-processor approach for the calculations. In this case, a simple application of OpenMP over 4 processors has yielded a threefold reduction in computation time for the particular system used by the author.
Mesh Quality Metrics
The three mesh quality metrics discussed here are the minimum/maximum included angle check, the aspect ratio check, and the warpage check.
The maximum angle check seeks to determine if any included angles (angles between two adjacent sides) of the quad element exceed a certain limit. For quad elements it is common practice to set this angle to 135°. A minimum angle check seeks to determine if any included angle is samller than a limit that is typically set to 45°.
The aspect ratio check seeks to determine the ratio of the length of largest element side to the length of the smallest element side. It is common practice to set the maximum allowable aspect ratio to 5.
The warpage check for quad elements seeks to assess the extent to which any one of the four nodes is not coplanar with the remaining three. The out-of-planeness is determined by dividing the quad element along a diagonal into two triangles, each of which defines a plane. The angle between the planes is measured. If this angle is greater than (for example) 15°, the quad element is deemed to have failed. The quad element is then divided into two different triangles by selecting the opposite diagonal as the dividing line. Once more, if this angle is greater than 15° then the quad element is deemed to have failed.
Sample Mesh
The sample mesh with poor quality elements that is used for this study is available for download here:
The mesh was generated using a crude mapped meshing technique and in all likelihood would never be used for a plane stress analysis. However, it does serve as a good test case for software development. The data within the file is stored in binary, little-endian format. The procedure for reading the file is quite straightforward:
1 Read two 4-byte integers - these should be 1 and 73344 (the number of nodes).
2 Read 73344 sets of (nodal) co-ordinates - three 8-byte doubles.
3 Read two 4-byte integers - these should be 4 and 72064 (the number of quads).
4 Read 72064 sets of (quad) vertices - four 4-byte integers per element.
The quad vertices are zero-indexed. Each index refers to one of the 73344 nodes.
When the mesh quality measurement is performed on the example mesh given in the above link, it can be seen that 15848 elements fail the minimum angle test with the worst angle being 22.48°. Furthermore, 15432 elements fail the maximum angle test with the worst angle being 157.52°. There are no aspect ratio failures and no warpage failures (which for the latter case is rather obvious owing to the fact that the mesh is contained in a single plane).
Source Code
The source code is available for download here:
The source code frequently uses objects of the V3 class, defined here in V3.h and implemented in V3.cpp. A V3 object is a 3d vector of type
double complete with a number of useful mathematical member functions such as addition, subtraction, scalar product, vetor product. These operations
are defined using operator overloading to provide a more intuitive read. Here are some typical exmaples of their use:
V3 u = { 0.5707, 0.5707, 0.5707 };
V3 v = { 1.0, 0.0, 0.0 };
V3 w;
w = u + v; // addition
double d = u * v; // dot product
w = u ^ v; // cross product
Note that, as is the case with a mathematical vector, a V3 object is used to contain a direction as well as a location.
The element quality check itself is performed by the function:
int ElementQualityCheck(int inodes[4], V3 xyz[], double vals[4]);
for which xyz is an array of V3 objects representing the locations of all the nodes in the mesh and
inodes is an array of 4 ints, each of which points to a node listed in xyz which together define
a quadrilateral element. On exit from the function, the array dvals of four doubles contains the maximum aspect ratio, the cosine of the maximum angle, the cosine of the minimum angle
and the cosine of the warpage angle respectively for those metrics for which the element has failed. If the element has passed
a metric, then the relevant entry in dvals will contain an 'idealized' value. This will assist in determining the
number of failed elements for a given metric.
dvals[0] = 1.0; // Idealized maximum aspect ratio.
dvals[1] = 0.0; // Idealized cosine of the maximum angle.
dvals[2] = 0.0; // Idealized cosine of the minimum angle.
dvals[3] = 1.0; // Idealized cosine of the warpage angle.
The function returns the number of quality metrics for which the element has failed (0 to 4).
A basic implementation is shown in the code snippet to the right in which the element check is wrapped in a loop.
The number of elements that have failed any of the quality metrics is summed in nfailed. The total number
of elements that have failed the minimum and maximum included angle check are assigned to nminang and nmaxang respectively.
The cosines of the maximum and minimum included angles encountered over all the elements are assigned to cminang and cmaxang respectively.
The code snippet assumes that the nodal locations and elemental definitions are stored in the following arrays:
V3 xyz[numNodes];
int inodes[numElements][4];
which will likely be dynamically allocated. The totoal number of elements to be analyzed is given by the int variable numElements.
Parallel Processing
For large meshes, the determination of quality for all the elements can be quite time consuming. Thankfully, the mesh analysis procedure described here uses a simple
for loop, with each iteration requiring roughly the same amount of work. This is a scenario which neatly lends itself to a multi-processor approach.
For the purpose of this example, some basic features of OpenMP are employed to divide the work evenly amongst a team of 4 threads (see source code to
the right). Work sharing is achieved using the loop construct #pragma omp parallel for.
The summation of the number of failed elements (nfailed) along with the summations of the number of elements that have
failed the minimum and maximum angle checks (
nminang and nmaxang respectively) are undertaken with the use of a
reduction clause that employs the + operator. Similarly,
the determination of the maximum and minimum included
angle that is encountered over all of the elements is achieved using a reduction clause that employs the min
and max operator. It should be noted that the min and max operators
were introduced in OpenMP 3.1, a standard which at the time of writing has not been adopted in every major C++ compiler. An alternative sample of source code (see right) is shown
in order to address such a situation.
Using parallel processing with 2 cores and 4 logical processors available from an Intel® Core™ i7-7500U Processor, the time required to check the quality of all elements in the sample mesh was reduced from .192 seconds to .071 seconds, which amounts to almost a threefold increase in performance.
double dvals[4];
int nfailed = 0;
int nminang = 0;
int nmaxang = 0;
double cminang = 0.0;
double cmaxang = 0.0;
for (int i = 0; i < numElements; ++i) {
if (ElementQualityCheck(inodes[i], xyz, dvals)) {
nfailed++;
cmaxang = min(cmaxang, dvals[1]);
cminang = max(cminang, dvals[2]);
nmaxang += (int)(dvals[1] < 0.0);
nminang += (int)(dvals[2] > 0.0);
}
}
omp_set_num_threads(4);
#pragma omp parallel for private(dvals) reduction(+:nfailed), \
reduction(+:nminang), reduction(+:nmaxang), \
reduction(min:cmaxang), reduction(max:cminang)
for(int i = 0; i < numElements; ++i) {
if (ElementQualityCheck(inodes[i], xyz, dvals)) {
nfailed++;
cmaxang = min(cmaxang, dvals[1]);
cminang = max(cminang, dvals[2]);
nmaxang += (int)(dvals[1] < 0.0);
nminang += (int)(dvals[2] > 0.0);
}
}
omp_set_num_threads(4);
double cmax[4] = { 0.0, 0.0, 0.0, 0.0 };
double cmin[4] = { 0.0, 0.0, 0.0, 0.0 };
#pragma omp parallel for private(dvals) reduction(+:nfailed), \
reduction(+:nminang), reduction(+:nmaxang)
for(int i = 0; i < numElements; ++i) {
if (ElementQualityCheck(inodes[i], xyz, dvals)) {
nfailed++;
cmax[omp_get_thread_num()] = min(cmax[omp_get_thread_num()], dvals[1]);
cmin[omp_get_thread_num()] = max(cmin[omp_get_thread_num()], dvals[2]);
nmaxang += (int)(dvals[1] < 0.0);
nminang += (int)(dvals[2] > 0.0);
}
}
cmaxang = cmax[0];
cminang = cmin[0];
for (int i=1; i < 4; ++i) {
cmaxang = min(cmax[i], cmaxang);
cminang = max(cmin[i], cminang);
}
Copyright © 2018 JGX Software Solutions LLC. All Rights Reserved.